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We provide a tight-binding model parametrization for black phosphorus (BP) with an arbitrary 
number of layers. The model is derived from partially self-consistent GWo approach, where the 
screened Coulomb interaction Wo is calculated within the random phase approximation on the basis 
of density functional theory. We thoroughly validate the model by performing a series of benchmark 
calculations, and determine the limits of its applicability. The application of the model to the 
calculations of electronic and optical properties of multilayer BP demonstrates good quantitative 
agreement with ab initio results in a wide energy range. We also show that the proposed model can 
be easily extended for the case of external fields, yielding the results consistent with those obtained 
from first principles. The model is expected to be suitable for a variety of realistic problems related 
to the electronic properties of multilayer BP including different kinds of disorder, external fields, 
and many-body effects. 

PACS numbers: 73.22.-f, 74.20.Pq, 71.10.Fd, 78.20.Bh 


I. INTRODUCTION 

A few-layer black phosphorus (BP) is a novel two- 
dimensional (2D) semiconductor with a number of re¬ 
markable properties such as strong anisotropy and pro¬ 
nounced thickness dependence of its electronic charac¬ 
teristics, which, along with high current on-off ratios and 
high carrier mobilities, make this material a promising 
candidate for diverse electronic and optical applications 
[IH3. Apart from the practical aspect, there is a grow¬ 
ing fundamental interest in BP ranging from attempts to 
provide insight into the origin of its band properties @ 
to more exotic and speculative aspects including super¬ 
conductivity and topologically nontrivial phases [T^ . 

From the theoretical perspective, one can distinguish 
between the two main approaches for studying electronic 
properties in material science. The first one is parameter- 
free first-principles calculations, commonly based on den¬ 
sity functional theory (DFT) and its many-body exten¬ 
sions (e.g., GW approximation). Although such meth¬ 
ods generally provide accurate results with respect to the 
ground state properties, their applicability to large sys¬ 
tems is very limited due to high computational cost and 
poor scalability. At the same time, realistic modeling 
in many cases requires large-scale simulations in order 
to, for example, describe finite-size effects, the presence 
of interfaces, or different kinds of disorder. In this re¬ 
spect, tight-binding (TB) Hamiltonian models act as an 
alternative approach to the electronic structure problem, 
providing a way to perform simulations with millions of 
atoms involved. Apart from being computationally more 
tractable, TB models also serve as a playground for ex¬ 
ploring rich many-body physics. 

Unlike graphene m , whose electronic properties in the 
low-energy limit are determined by a simple TB Hamil¬ 
tonian, involving only one nonequivalent parameter (in¬ 
tersite hopping integral, t), a reliable theoretical descrip¬ 


tion of a single-layer BP (known as phosphorene) is con¬ 
siderably more challenging. A number of low-energy 
electronic properties of pristine phosphorene can be effi¬ 
ciently described in terms of the (2x2) k ■ p Hamilto¬ 
nian [a, EMI with parameters determined to reproduce 
first-principles calculations. Being determined in recip¬ 
rocal space and designed to describe the valence band 
(VB) and conduction band (CB) edges only, the k • p 
Hamiltonians are not well suited for studying real-space 
problems. Moreover, such models basically rely upon the 
effective mass approximation, whose applicability is not 
well justified for BP even in the low-energy range due to 
the presence of flat bands. Last but not least, although 
the extension of the k • p model appears straightforward 
to the multilayer case 011113 , it becomes dependent 
on thickness-dependent parameters, which are a priori 
not known. 

Early attempts to provide a real-space model to the 
electronic structure of BP were based on molecular or¬ 
bital theory [Tl |. whose simplified nature and complex 
orbital character of BP prevent a quantitatively accu¬ 
rate description [I^. Recently, two of us have proposed 
a more rigorous real-space model for single- and double¬ 
layer BP, which was constructed by downfolding the full 
Go Wo Hamiltonian to the minimal (one interaction site 
per phosphorus atom) low-energy effective Hamiltonian 
[20|. The latter involves two main parameters of un¬ 
like signs corresponding to two nearest-neighbor hopping 
integrals, and a number of less-relevant long-range pa¬ 
rameters needed to accurately reproduce the quasiparti¬ 
cle VB and CB edges of monolayer and bilayer BP. The 
model has been successfully applied in a number of stud¬ 
ies including those related to phosphorene nanoribbons 
min, electric and magnetic fields EH, dif¬ 

ferent kinds of disorder [EJl, and realistic modeling of 
field-effect electronic devices [13 . However, the appli¬ 
cability of that model is limited to single- and bilayer 
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BP, whereas thicker (experimentally available) samples 
cannot be considered. 

In this paper, we report on a revision of the above men¬ 
tioned model [^. Particularly, we focus on its modifi¬ 
cation to describe BP samples with arbitrary thickness, 
ranging from monolayer to bulk. We also improve the 
quantitative validity of the model, which allows us to 
achieve consistency with experimental results in the bulk 
limit. The proposed model is derived on the basis of accu¬ 
rate first-principles calculations within the partially self- 
consistent GWq approximation and systematically vali¬ 
dated by performing a series of benchmark tests. The 
model is suitable for studying large-scale problems and 
applicable in a wide energy range. As a study case, we 
examine the energy gap dependence on the number of 
layers and also consider the influence of a perpendicular 
electric field onto the electronic structure of BP. Particu¬ 
larly, we study the role of BP thickness in the transition 
from a normal to topological insulator driven by external 
electric field recently predicted for a few-layer BP 0 - 

The rest of the paper is organized as follows. In Sec. II, 
we start with an overview of previous first-principles 
studies of the electronic structure of BP (Sec. II A), pro¬ 
vide calculation details (Sec. IIB), and present the results 
of the GWo calculations, accompanied by the analysis of 
the quasiparticle band structure of a few-layer and bulk 
BP (Sec. IIC). In Sec. Ill, we propose the TB model, 
describe the parametrization procedure (Sec. Ill A), and 
perform a series of calculations in order to assess its per¬ 
formance (Secs. IIIB and IIIC). In Sec. IV, we extend 
the model by adding an electric field and apply it to mul¬ 
tilayer BP. In Sec. V, we briefly summarize our results. 


II. ELECTRONIC STRUCTURE OF A 
FEW-LAYER BP FROM FIRST PRINCIPLES 

A. Overview of previous studies 


After a few-layer BP became available experimentally, 
a considerable number of theoretical studies of its band 
structure have been reported. The calculations showed 
that commonly used DPT in conjunction with local and 
semilocal exchange-correlation approximations does not 
describe semiconducting properties of bulk BP correctly. 
Contrary to experimental observations, yielding a nar¬ 
row gap of 0.31-0.35 eV [26l - l^ for bulk BP, the local 
density approximation or generalized gradient approxi¬ 
mation (GGA) predict signihcantly smaller or even zero 
values, depending on a particular computational scheme 
and lattice parameters |^l29l - l34 | . The utilization of hy¬ 
brid functionals [such as Heyd-Scuseria-Ernzerhof (HSE) 
HH], ncorporating a nonlocal contribution to the ex¬ 


act exchange, has been shown to partially solve the band 
gap problem of bulk BP [29l - [3lll^l34 1. However, the per¬ 
formance of such methods depends strongly on a number 
of empirical parameters determining, for example, the 
screening range and fraction of the exact exchange con¬ 


tribution, which are generally system specific and cannot 
be systematically determined. This ambiguity results in 
a broad variation of band gaps in a few-layer and bulk 
BP (see Table U for an overview). 

More consistent results with respect to the band prop¬ 
erties can be obtained using the GW a ppr oximation |37l| , 
which has been applied to BP in Refs. |20ll3^l33l^ . The 
authors of Ref. adopt a non-self-consistent GqWq 
scheme, where the screened Goulomb interaction Wo is 
calculated within the general plasmon pole model 
and report a band gap of 0.3 eV for bulk BP, which is 
within the range of available experimental data. How¬ 
ever, the use of a more reliable random phase approx¬ 
imation (RPA) [ 4 ^ within the GqWq scheme yields a 
smaller value of 0.1 eV [ 2 ^. More accurate band gap val¬ 
ues are supposed to be obtained within the RPA in terms 
of a partially self-consistent GWo scheme. Such calcula¬ 
tions have been recently performed in Ref. , where the 
evaluation of Wq was based on hybrid functionals and re¬ 
sulted in significantly higher band gap values for bulk BP 
(0.58 eV) compared to the experimental ones. Therefore, 
the hybrid functionals do not seem to be an optimal start¬ 
ing point for GW calculations of BP. Physically, this can 
be attributed to excessively contracted wave functions, 
which suppress the screening of the Goulomb repulsion 
and eventually leads to the band gap overestimation. The 
closest results to experiment are obtained by means of the 
GWo approach with Wo calculated on top of the GGA 
wave functions within the RPA (denoted as GITo@GGA 
thereafter) [ 3 ^, which yields the gaps of 0.43 and 1.94 
eV for bulk and monolayer BP, respectively. The latter 
value is also consistent with recent scanning tunneling 
spectroscopy measurements of the gap in the spectrum 
of surface states of cleaved BP (2.05 eV) (3^ . 


B. Calculation details 

Here, we first apply the GITo@GGA scheme to calcu¬ 
late the quasiparticle electronic band structures for n- 
layer (n = 1-3) and bulk BP, which provide reference 
data for the subsequent TB model parametrization. The 
calculations were performed within the projected aug¬ 
mented wave formalism [4ll| as implemented in the Vi¬ 
enna ab-initio simulation package (vASP) (^|4^ . The 
Green’s functions (G) were first calculated by using the 
Kohn-Sham eigenvalues and eigenstates and then iter¬ 
ated four times, which proved to be sufficient to achieve 
numerical convergence |44l | . The screened Goulomb inter¬ 
action (Wo) is calculated on the basis of the frequency- 
dependent dielectric function, Wn = Cq which, in turn, 
is computed at the RPA level [4^ as cq = I — vxo, where 
V is the bare Goulomb interaction and xo is the inde¬ 
pendent particle polarizability. The latter is evaluated 
by using the DFT-GGA [i^ eigenvalues and eigenstates 
in the spectral representation. To this end, a numerical 
integration along the frequency axis containing 70 grid 
points is performed. In the calculation of the quasiparti- 
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FIG. 1. (Color online) Orbitally resolved densities of states (DOS) calculated for a few-layer (n = 1-3) and bulk (n = oo) BP by 
projecting the GWo Hamiltonian onto the atom-centered Wannier orbitals, corresponding to s and pi (i = x,y,z) symmetries. 
The total DOS is shown in the inset within the energy range of —8 to 8 eV relative to the band gap center indicated by the 
vertical dashed line. 


cle energies, both diagonal and off-diagonal elements of 
the self-energy matrix S = iGW were included. The to¬ 
tal energy in the DFT part was converged to within 10“® 
eV. In all calculations, we use an energy cutoff of 250 eV 
for the plane-wave expansion of the wave functions. The 
number of unoccupied states in GW calculations were 
set to 90 per atom. In most cases, a k-point mesh of 
(10 X 12 X 1) and (10 x 12 x 4) was used for the Brillouin 
zone sampling of a few-layer and bulk BP, respectively. 
To examine the fine structure of the electronic spectrum 
of monolayer BP, a denser mesh was considered. To ob¬ 
tain smooth band structures, densities of states and op¬ 
tical conductivities, we use an interpolation procedure 
by makin g u se of the maximally localized Wannier func¬ 
tions which are constructed by projecting the 

GWq Hamiltonian onto the entire manifold of the 3s and 
3p states of phosphorus. For all the structures, we adopt 
experimental crystal structures of bulk BP and in¬ 
troduce a vacuum layer of ~15A in order to minimize 
spurious effects due to the periodic boundary conditions 
in slab calculations. The chosen set of parameters ensures 
that the quasiparticle gaps are accurate to within a few 
hundredths of eV. Although some variations in structural 
parameters have been reported between monolayer and 
bulk BP [mi, we intentionally do not consider such 
effects in our work due to the following reasons: (i) to 
minimize the complexity of the TB model for multilayer 
BP associated with atomic degrees of freedom, and (ii) 
to avoid ambiguity in the determination of structural pa¬ 
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FIG. 2. (Golor online) Left: Fine structure of the VB of mono¬ 
layer BP calculated along the V-Y direction in the vicinity of 
the F point by means of the DFT-GGA and GWa approaches. 
Points correspond to the original calculations at a (24 x 32 x 1) 
k-point mesh, whereas lines represent Wannier-interpolated 
bands. Right: Wannier-interpolated DFT-GGA and GWq 
densities of states calculated for the same energy range. Zero 
energy corresponds to the position of the VB at the F point. 


rameters for a few-layer BP at the DFT level arising from 
the variety of different exchange-correlation functionals. 


C. Quasiparticle electronic properties from 
partially self-consistent GWq approximation 

In Fig. [U we show the densities of states (DOS) cal¬ 
culated for a few-layer (n = 1-3) and bulk BP within 
the GlFo@GGA scheme. One can see that the calculated 
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value of a band gap of bulk BP is 0.35 eV, which is within 
the bounds of experimental variability (0.31-0.35 eV [2^ 
Hi). Such an agreement justifies the computational ap¬ 
proach employed and allows us to expect accurate results 
for a few-layer BP. Qualitatively, the GWq results for a 
few-layer BP shown in Fig.[T]are similar to those reported 
previously (20ll30ll33 - l^ | . In cases of a few-layer BP, DOS 
exhibits a step like behavior, which is typical for systems 
with reduced dimensionality. For the following, it is also 
worth mentioning that in all cases considered, the major 
contribution to the states close to the band gap comes 
from the Pz states of phosphorus, whereas Py states have 
zero contribution at the VB and CB edges. 

In the case of monolayer BP, the fine structure of the 
electronic states near the edge of the VB requires special 
attention. As has been noticed in previous DFT stud¬ 
ies [l^, the VB maximum is slightly shifted from the F 
point in the F-V direction, which apparently results in 
an indirect gap in monolayer BP. The deviation of the 
VB maximum from the zone center might result in non¬ 
trivial physical properties of BP such as superconducting 
and ferromagnetic instabilities [s^ due to the appear¬ 
ance of the van Hove singularity close to the VB edge. 
At the level of the k • p perturbation theory, a transition 
from a direct to an indirect band gap in monolayer BP 
is governed by the magnitudes of the matrix elements 
of the momentum operatorjCorresponding to transitions 


between the VB and CB 


At the same time, well- 


known inaccuracies of DFT with respect to the VB and 
CB positions cannot support the prediction of an indirect 
gap in monolayer BP. Therefore, it appears appropriate 
to examine the fine structure of the monolayer VB at 
the more accurate GWq level. To this end, we perform 
a comparison between the electronic structures of mono- 
layer BP calculated within the DFT-GGA and GWq ap¬ 
proaches by using a dense (24 x 32 x 1) k-point mesh. 
The results are shown in Fig. [21 We do reproduce the 
previously reported shift of the VB maximum from the 
F point as well as the van Hove singularity in DOS calcu¬ 
lated at the DFT-GGA level. However, the GWq results 
show no indications of such a behavior and support for a 
direct band gap in monolayer BP. 


III. TIGHT-BINDING MODEL FOR 
MULTILAYER BP AND ITS VALIDATION 

A. Parametrization procedure 

The effective TB model considered in this work is given 
by the effective four-band Hamiltonian, describing one 
electron per lattice site, 

^ ( 1 ) 

where i and j run over the lattice sites, (t^) is the 
intralayer (interlayer) hopping parameter between the i 



FIG. 3. (Color online) Schematic representation of the hop¬ 
ping parameters for the TB model [Eq. (|T])] parametrized in 
this work for multilayer BP. The corresponding values are 
given in Table [III 


and j sites, and c] (cj) is the creation (annihilation) oper¬ 
ator of electrons at site i (j). We note that in contrast to 
the model used in our previous works [^1^ , the Hamil¬ 
tonian given by Eq. m does not contain on-site terms, 
meaning that electrons at all sites have equivalent ener¬ 
gies even for multilayer BP. 

To parametrize the model given by Eq. (HD , we use a 
procedure similar to Ref. [2^, which is as follows. We first 
map the entire manifold of valence and conduction states 
of phosphorus monolayer onto the subspace of effective 
Pz-like orbitals (four orbitals per unit cell) being relevant 
for the low-energy part of the VB and CB. To this end, 
we hrst use the original Bloch states |^nk) obtained from 
the GWq calculations and construct a new subspace of 
Bloch-like states IV'nk), 

p 

|^„k) = ^ UlMm,.), (2) 

m—l 

where P is the total number of states included into the 
GWo calculations and is a rectangular matrix ob¬ 
tained by projecting the \pz) states onto the Bloch states 
l^/^nk) and using the disentanglement procedure proposed 
in Ref. [5l|. Having obtained I'i/'nk), we construct an ef¬ 
fective (4 X 4) Hamiltonian in reciprocal space iJm„(k), 
which is achieved by performing a unitary transforma¬ 
tion of the original GWq Hamiltonian in the Bloch 
subspace. The resulting reciprocal-space Hamiltonian 

= (^mk|i?1^nk) (3) 

is then transformed into the real space, = 

The resulting real-space Hamiltonian 
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is determined in the basis of Wannier functions 
1 '^^) = Sk e”**' ^iV'nk), corresponding to the p^-like or¬ 
bitals. 

Despite low dimensionality of its matrix elements 
(hopping parameters) decay slowly with distance, result¬ 
ing in a large number of small parameters. In order to 
make the resulting model more tractable, we ignore the 
parameters beyond the cutoff radius of ~5.5 A, which are 
typically smaller than 0.01 eV. To restore the quality of 
the truncated Hamiltonian, we reoptimize the remaining 
parameters in such a way that they provide an accurate 
description of the band structure in the low-energy re¬ 
gion. To this end, we minimize the following least squares 
functional, F{{U}) = 

where {ti} are hopping parameters and (e™) 

an eigenvalue of the corresponding (GWo or TB model) 
Hamiltonian n and k are the band index and 

momentum vector, respectively, which run over the rele¬ 
vant region in the vicinity of the band gap. In the case 
of monolayer, this region involves the valence and con¬ 
duction bands only. To parametrize the TB Hamiltonian 
for bilayer, we adopt a similar strategy. In this case, we 
introduce interlayer hopping parameters, while the in¬ 
tralayer parameters remain fixed. Also, we take into ac¬ 
count the splitting of the valence and conduction bands 
upon the optimization of the hoppings, which is crucially 
important for the applicability of the model to multilayer 
BP. The obtained set of intralayer ({tf}) and interlayer 
hoppings are then applied without any corrections 
to BP with a larger number of layers. 


B. Electronic structure 

The resulting hopping parameters are schematically 
shown in Fig.[3]and listed in Table HIl Overall, our model 
involves ten intralayer and four interlayer hoppings. As 
has been previously noticed [^1^ , the main features of 
the band structure of monolayer BP can be qualitatively 
described by only two largest hopping parameters (tf and 
t|). The band gap at the P point is determined in this 
case by a simple expression, « 2\t\ \ — 4|tf |. For 

bilayer BP, the degeneracy of the VB and CB is lifted if 
a nearest-neighbor interlayer hopping (tj-) is introduced. 
This results in a reduction of the band gap, given now 

by£’f^(r) - - 4|tf I - 2\ti\. In order 

to quantitatively reproduce the quasiparticle spectrum of 
BP including accurate k dependence of the VB and CB 
as well as their splitting in the multilayer case, a larger 
number of hopping parameters is required. 

In Fig. m we show the band structures calculated 
within the derived TB model in comparison with the 
full bands obtained from GWq calculations. One can see 
that the TB model accurately describes the results of 
GWq calculations in the low-energy region not only for 
monolayer and bilayer BP, but also for trilayer and bulk 


structures. Since the band properties of trilayer and bulk 
BP have not been used as a reference during the model 
parametrization, it is natural to expect the applicability 
of the presented model to BP with an arbitrary number 
of layers. 

To explicitly demonstrate that the obtained TB Hamil¬ 
tonian is represented in a physically meaningful orbital 
subspace corresponding to the p^-like states, we consider 
the case of monolayer BP, for which we project the full 
GWq band structure onto the states [see Fig.jSfa)] and 
compare it with the model bands [Fig. [SJb)]. From the 
projected GIFo bands shown in Fig.jSJa) one can clearly 
recognize four distinct bands having predominantly Pz 
symmetry, whose contribution is shown by color. By 
comparing those with Fig. EDb) it becomes evident that 
the model provides an effective representation of the Pz- 
like states. As can be inferred from Fig. [5l[a) and will 
be shown below, the states of the other symmetries do 
not contribute to direct interband transitions within an 
energy range of up to several eV, which basically deter¬ 
mines the limits of the applicability of the presented TB 
model. 

Having obtained a TB model applicable for multilayer 
BP, it is instructive to analyze the the band gap de¬ 
pendence on the number of layers. In Fig. [51 we show 
the corresponding dependence calculated within the TB 
model, which can be accurately htted by the expression 
= Aexp{—nB)/n'^ + D with parameters A, B, G, 
and D given in the inset of Fig. [51 One can see that 
along with a power law decay, being important at small 
n, there is a pronounced exponential decay, becoming 
dominant at large n. Our result is thus different from 
the previously proposed power law expected from a sim¬ 
ple quantum confinement picture [s^. 

C. Optical properties 

To further validate our model, we calculate the 
frequency-dependent optical conductivity Uajsiuj) calcu¬ 
lated for the undoped case by means of the GWq ap¬ 
proach and TB model for n = 1-3 layer and bulk BP. 
Within the GWo, we evaluate aa 0 {oj) through the Bril- 
louin zone integration using the following form of the 
Kubo-Greenwood formula in the independent-particle 
approximation 

(J /mk - /nk {n'k\ja\rn'k){rnk\jf 3 \nk) 

, ^mk £^nk {hlO 4- ip) 

k mn 

( 4 ) 

where is the unit cell area, Nk is the number of k 
points used for the Brillouin zone sampling, |mk) is the 
Wannier-interpolated Bloch state [i^, corresponding to 
the mth eigenvalue Smk of the GWq Hamiltonian , 

/nk = exp(/3£„k + 1)”^ is the Fermi-Dirac occupation 
factor involving the inverse temperature /3, ja is the a 
component of the current operator, and ry is a smearing 
parameter. The Brillouin zone was sampled by ^10^ and 
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FIG. 4. (Color online) Electronic band structures calculated for a few-layer (n = 1-3) and bulk BP by using the GWq 
approximation (light lines) and within the TB model given by Eq. ([T]) with parameters listed in Table HIl (dark lines). Zero 
energy corresponds to the center of the gap. High-symmetry points of the Brillouin zone are shown in the insets. 



FIG. 5. (Color online) Band structure and density of states 
(DOS) calculated for monolayer BP by using (a) GWo ap¬ 
proach and (b) TB model of this work. In (a), the DOS are 
projected onto the Pz states, whereas in the band structure 
their contribution is shown by color. High-symmetry points 
of the Brillouin zone are shown in the insets. 


FIG. 6. (Color online) Layer dependence of a band gap in 
BP calculated by using the GWo approximation, TB model 
presented in this work and by an empirical expression = 
Hexp(—nB)/n‘^ + D with parameters A, B, G, and D fitted 
to the TB model. 


10® k points for 2D (a few layer) and 3D (bulk) calcu¬ 
lations, respectively. To demonstrate the advantage of 
the derived TB model for studying realistic samples, we 
apply the TB Hamiltonian [Eq. ([I])] to calculate aapioj) 
for a few-layer (n = 1-3) and many-layer (n=100) BP in 
real space. To this end, we use the tight-binding prop¬ 


agation method in which is calculated 

conceptually similar to Eq. o but by considering ex¬ 
plicit evolution of the current operator in time [e.g., see 
Eq. (30) of Ref. [ 5 ^] instead of diagonalization of large 
matrices. The sample size was taken to contain ^10^ 
atoms in each case considered with periodic boundary 
conditions applied in lateral {xy) directions. In both 
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methods, we restrict ourselves to the diagonal compo¬ 
nents of ao,p{uj) only. We stress that Uapiuj) is calculated 
within a single-particle approximation, meaning that the 
excitonic effects are neglected. Such effects are proven to 
be relevant for a reliable description of the optical spec¬ 
tra of monolayer and a few-layer BP, but they become 
insignificant in the bulk limit (32|. 

The results of our calculations are shown in Fig. [T] 
In line with previous studies [mi, we observe strong 
anisotropy between the conductivities in different direc¬ 
tions and well-pronounced peaks along the armchair di¬ 
rection of a few-layer BP, associated with the discrete 
character of the band structure close to the VB and CB 
edges. As can be seen from Fig.[7l the optical conductiv¬ 
ities obtained with the use of the TB model are in very 
good agreement with the results of GWq calculations in a 
wide frequency range up to 2.0 eV. The agreement in the 
range of 2.0-2.5 eV can be considered as satisfactory but 
it is becoming worse for structures with a large number of 
layers. At larger frequencies (w > 2.5 eV), the TB model 
still shows reasonable agreement with the GWq results for 
a few-layer BP, but becomes apparently inapplicable to 
many-layer systems (including bulk), which is due to the 
decreased gap, allowing for transitions between the states 
not included in the construction of the TB model. De¬ 
spite being relatively close to the band gap, those states 
do not contribute to the optical conductivity at the lower 
frequencies since the expression for aap{uj) [Eq. (|4])] in¬ 
volves only direct transitions between the VB and CB. 
We note that for transport and optical properties involv¬ 
ing indirect transitions between the VB and CB (e.g., in 
scattering processes) a reliable frequency range for the 
TB model will be more limited and determined entirely 
by the consistency between the quasiparticle and model 
bands shown in Fig. 2] 

IV. EFFECT OF ELECTRIC FIELD ON THE 
BAND STRUCTURE OF MULTILAYER BP 

We now consider an extension of our model to the case 
of an electric field Ez perpendicular to the surface. For 
simplicity, we restrict ourselves to bilayer BP, for which 
the extended Hamiltonian reads 

H = +eEzZ, (5) 

where the first term in the right-hand side corresponds to 
the unperturbed Hamiltonian for bilayer given by Eq. (HD 
and the second term plays the role of a layer-dependent 
on-site potential. We note that in what follows, we con¬ 
sider an unscreened electric field only, that is, we neglect 
explicit treatment of polarization and local-field effects. 
In other words, Ez can be regarded as a local electric 
field assumed to be constant inside the sample. Ez can 
be related to real external electric field A™* upon tak¬ 
ing into account thickness-dependent transverse dielec¬ 
tric permittivity Szid) and finite-size effects. In a first 
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FIG. 7. (Color online) Optical conductivities for a few-layer 
(n = 1-3) and bulk BP calculated along the armchair 
and zigzag (oyy) directions using the Kubo formula [Eq. ([4|l] 
on the basis of the GWo approach and TB model presented 
in this work. (Jz:x:(yy) are given per layer in terms of the uni¬ 
versal optical conductivity of graphene (no = jAh). Within 
the TB model, bulk BP is approximated by 2D BP with a 
large number of layers (n =100). In all cases, we set the 
temperature to 300 K. 


approximation, one can take = EzEz, where £z is the 
transverse dielectric permittivity of bulk BP (e^ 8.3 

H). 

In Fig. [SI we show the low-energy part of the band 
structure calculated for three representative electric 
fields. In the presence of an electric field, the electronic 
bands shift due to the difference of the interlayer poten¬ 
tial which is a manifestation of the Stark effect. From 
Fig. ISI one can see that the VB and CB shift in differ¬ 
ent directions toward the band gap center. This causes 
a decrease of the band gap with increasing field, which 
reaches zero at Ez = 341 mV/A. At higher field the band 
inversion is observed, as can be seen from Fig. |Sl[f). Our 
results obtained using the TB model are thus consistent 
with previous DFT calculations for a few-layer BP and 
phosphorene nanoribbons [miii-iiii. 

It is interesting to note the existence of a Dirac-like 
linear dispersion along the armchair direction (A-F-A) 
at the critical electric field, E'^ [Fig. |5Kc)], which appears 
around the F point. A qualitatively different situation is 
observed in the zigzag direction (Y-F-V), where the dis¬ 
persion turns out to be quadratic [Fig. |51[d)]. At higher 
fields {Ez > Ef) the Dirac point disappears in the arm¬ 
chair direction, whereas two band crossings appear along 
the zigzag direction as a result of the band inversion. 

Finally, we calculate the evolution of the critical bias 
potential, = eE‘^d with the number of BP layers n, 
which is applied between the top and bottom planes of an 
n-layer sample separated by the distance d. In Fig.jHl the 
corresponding dependence is shown. Since is related 
to the original band gap at zero field, it is natural to 
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FIG. 8. (Color online) Band structures of bilayer BP calcu¬ 
lated in the vicinity of the F point for different magnitudes of 
the electric field = E% + AEz, where E^ =341 mV/A is 
a critical field at which the band gap closes, and AEz takes 
the values of —2, 0, and +2 mV/A. Top (bottom) panels cor¬ 
respond to the bands calculated along the X-T-X {Y-T-Y) 
directions. Valence and conduction bands are indicated by 
blue and orange, respectively. Zero energy corresponds to the 
center of the gap at the F point. 



FIG. 9. (Color online) Critical bias potential AV^(n) required 
to close the band gap in n-layer BP plotted as a function of 
the number of layers. The bias potential is applied within 
the TB model of this work [see Eq. ®] between the top and 
bottom planes of the corresponding BP samples neglecting the 
screening effects. Points correspond to the TB calculations, 
whereas lines correspond to the fitting via AV'^(n) = AjvP Y 
D. 


expect that the same form of the functional dependence 
as in Fig.[5](i.e., a power law with exponential cutoff) can 
be used to parametrize AV^{n). We find, however, that 


in the present case the prefactor B in the argument of 
the exponential is significantly smaller {B < 0.01). This 
allows us to fit the critical bias potential as AV^(n) = 
A/nP + D, where A, C, D are fitting parameters given in 
the inset of Fig. |9l We conclude, therefore, that AV^{n) 
(Fig. [S]) exhibit a significantly weaker dependence on the 

number of BP layers than the band gap, (Fig- El)- 
V. CONCLUSIONS 

We have proposed an effective TB model for multi¬ 
layer BP with arbitrary thickness, which is parametrized 
on the basis of partially self-consistent GWq approxima¬ 
tion. The model shows good performance with respect 
to static band properties as well as transport character¬ 
istics of multilayer BP compared to the GWq results. 
In contrast to previously proposed k • p Hamiltonians 
for BP, our model (i) directly applicable in real space; 

(ii) goes beyond the effective mass approximation; and 

(iii) accurately reproduces low-energy electronic proper¬ 
ties of multilayer BP without the need for additional 
scaling parameters. On the other hand, the proposed 
model is substantially less computationally demanding 
than any first-principles calculations, which makes cal¬ 
culations with millions of atoms possible. This allows 
us to expect its suitability for use in investigations of 
a wide range of phenomena, particularly in large-scale 
simulations of realistic BP (e.g., with disorder or in the 
presence of external fields) and as a starting point for 
studying many-body effects in BP. As an example of the 
model extension, we considered the case of an electric 
field applied to multilayer BP, which allowed us to deter¬ 
mine the thickness dependence of the critical bias poten¬ 
tial required to reach the regime of the band inversion 
previously predicted by first-principles calculations. We 
also found that the critical bias potential decays signif¬ 
icantly slowly with the number of BP layers than does 
the corresponding band gap. 
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TABLE 1. Band gaps (in eV) for monolayer (n = 1), multilayer (n = 2, 3), and bulk BP (n = c») calculated at different levels 
of theory. In the notation of different methods, Go and Wo imply that the Green’s function and screened Coulomb repulsion in 
the GW approach are calculated non-self-consistently on the basis of wave functions derived from density functional (GGA) or 
hybrid functional (HSE) calculations, whereas G means a self-consistent calculation of the Green’s function. Wp and Wo denote 
that the screened Coulomb interaction is calculated by using the general plasmon pole model [S^ and RPA [^, respectively. 




GIVo@GGA^ 

TB ModeP 

GITo@GGA^ 

GITo@HSE^ 

GoWo@GGA^ GoIPo@GGA^ 

HSE^ 

GGA^ 

Expt. 

n 

= 1 

1.85 

1.84 

1.94 

2.41 

1.60 

2.00 

1.00-1.91 

0.80-o”l 

2. 05^1 

n 

= 2 

1.16 

1.15 

~1.65 

1.66 

1.01 

~1.30 

1.01-1.23 

0.45-0.60 

— 

n 

= 3 

0.84 

0.85 

~1.35 

1.20 

0.68 

~1.05 

0.73-0.98 

0.20-0.40 

— 

n - 

= oo 

0.35 

0.40 

0.43 

0.58 

0.10 

0.30 

0.18-0.39 

0.00-0.15 

0.31-0.35^ 


® This work. 
Reference 
Reference 
Reference 
® Reference 
^ References 



^ References _ _ 

This value corresponds to a gap in the spectrum of surface states of bulk BP (Ref. l38l l. 
‘ References [2^428 


TABLE 11. Intralayer (t'l) and interlayer (t^) hopping parameters (in eV) obtained in terms of the TB Hamiltonian [Eq. fT])] 
for multilayer BP. d and Nc denote the distances between the corresponding interacting lattice sites and coordination numbers 
for the given distance, respectively. The hoppings are schematically shown in Fig. [S] 



Intralayer 



Intralayer 



Interlayer 


No. 

t" (eV) 

d(A) 

Np 

No. 

t" (eV) 

d(A) 

Na 

No. 

t- (eV) 

d(A) 

iVc 

1 

-1.486 

2.22 

2 

6 

0.186 

4.23 

1 

1 

0.524 

3.60 

2 

2 

3.729 

2.24 

1 

7 

-0.063 

4.37 

2 

2 

0.180 

3.81 

2 

3 

-0.252 

3.31 

2 

8 

0.101 

5.18 

2 

3 

-0.123 

5.05 

4 

4 

-0.071 

3.34 

2 

9 

-0.042 

5.37 

2 

4 

-0.168 

5.08 

2 

5 

-0.019 

3.47 

4 

10 

0.073 

5.49 

4 

5 

0.000 

5.44 
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